python3 scripts/convertToGFA.py {input} {output} {params.k} 2> {log}

fastp -i {input} -o {output.fixedFasta} -p -j {output.reportjson} -h {output.reporthtml}

import matplotlib.pyplot as plt
import numpy as np
import scipy.signal
import json

kmer_counts = json.load(open(snakemake.input['kmerProfile'],'r'))
uniqueKmersPerK = json.load(open(snakemake.input['uniqueKmers'],'r'))
k= snakemake.params['k']

plt.rcParams['figure.figsize'] = [16, 24]

fig,axs = plt.subplots(2)

kmer_counts_list = [kmer_counts[kmer] for kmer in  kmer_counts]
b = np.amax(kmer_counts_list)

hist,_,_ = axs[0].hist(kmer_counts_list, bins=b)
axs[0].set_title('raw k-mer histogram')

local_mins = scipy.signal.find_peaks_cwt(-hist,np.arange(30,70))
threshold = local_mins[0]


# Create a list of k-mers that occur more frequent than the cutoff-point
filteredKmers = [kmer for kmer in kmer_counts if kmer_counts[kmer] > threshold]
kmer_counts_list_filtered = [kmer_counts[kmer] for kmer in filteredKmers]
b = np.arange(np.amax(kmer_counts_list_filtered) + 1)
hist_filtered, _, _ = axs[1].hist(kmer_counts_list_filtered, bins=b)

# Compare to expected value
axs[1].set_title(
    'error-threshold: {} | k: {} | unique k-mers: {}/{} expected'.format(threshold,k,len(filteredKmers),uniqueKmersPerK[k])
,fontsize=9)


    
fig.text(0.5,0.04,'k-mer frequency',ha='center')
fig.text(0.04,0.5,'k-mer count',va='center',rotation='vertical')
plt.savefig(snakemake.output[0])

from Bio import SeqIO
from shared import *
import matplotlib.pyplot as plt
import json

x = []
y = []

uniqueKmersPerK = {}

# Read reference
reference = SeqIO.read(snakemake.input[0],'fasta')

ks = snakemake.params['krange']

for k in ks:

    # Create an empty dictionary to score kmer counts (how often each kmer is observed)
    kmers_reference = {}

    parseKmers(kmers_reference,reference,k)
    
    
    ambiguousKmers = sum(1 for kmer in kmers_reference if kmers_reference[kmer] > 1)
    totalKmers = sum(1 for kmer in kmers_reference)
    ratio = ambiguousKmers/totalKmers
    
    uniqueKmersPerK[k] = totalKmers-ambiguousKmers
    

    
    duplicates = sum(1 for kmer in kmers_reference if kmers_reference[kmer] == 2)
    
    print(k,totalKmers,duplicates)
    
    x.append(k)
    y.append(ratio)
  
# plot tobigram
  
plt.rcParams['figure.figsize'] = [10, 5]
plt.plot(x,y)
plt.xlabel('k')
plt.ylabel('fraction of ambiguous k-mers')

plt.savefig(snakemake.output['tobigram'])
with open(snakemake.output['uniqueKmers'],'w') as outfile:
	json.dump(uniqueKmersPerK,outfile)
    

bcalm -in {input} -kmer-size {params.k} -out {params.outprefix} -abundance-min 50 2> {log}

perl scripts/fasta_to_fastq.pl {input} > {output}

from Bio import SeqIO
from shared import *
import json

scReads = SeqIO.parse(snakemake.input[0],'fasta')

k = snakemake.params['k']

kmers = {}

for read in scReads:
	parseKmers(kmers,read.seq,k)

with open(snakemake.output[0],'w') as outfile:
	json.dump(kmers,outfile)

import pysam
from Bio import Seq,SeqIO,SeqRecord

alignment = pysam.AlignmentFile(snakemake.input['alignment'],'rb')
records = []

for segment in alignment.fetch():
	seq = Seq.Seq(segment.query_alignment_sequence)
	rec = SeqRecord.SeqRecord(seq, segment.qname, "", "")
	records.append(rec)

with open(snakemake.output[0],'w') as outfile:
	SeqIO.write(records,outfile,'fasta')

samtools index {input}